REM------- Q2ma ----- a Rational function(2nd order)-------------------------
REM
REM  Mandelbrot set of
REM             Q(z)=z/(z^2-az+a),   z=x+i*y,  a=ar+i*ai
REM
REM  2つのcritical points c1=√a  and c2=-√a  を初期値 Z0にしてのiteration
REM 
REM             Zn+1=Q(Zn),  n=0,1,2,....
REM               O+(Z0): Forward Orbit of Z0
REM 
REM -----------------------------------------------------------------
REM        カラー      ｜  O+(-√a）収束    ｜  O+(-√a）収束しない
REM ----------------------------------------------------------------- 
REM  O+(√a) 収束      ｜    co(1,1)黒 1    ｜    co(1,2)赤 4
REM  O+(√a) 収束しない｜    co(2,1)緑 3    ｜    co(2,2)黄 6
REM ----------------------------------------------------------------
REM        黒以外のカラーは変更可
REM ＊簡便法 
REM   O+(ci) が収束しない  <＝>  ∃n(>10); |Zj| >chao　(j>n)　となった時
REM (注) 10<= n <tm だが, n=10〜100位か（100のときはtm=400以上が良い）
REM       chao=3.5〜4 位でよい                                 
REM         
REM     iex(k), k=1,2  ... 分母ほとんど0の回数
REM     （分母ぴったり0になると，計算stopする => xl, xu を選び直せ）           
REM     tex(k)... 反復n回以下で絶対値 chao以上の回数
REM---------------------------------------------- Joe 5.31 '19 ---------------         
DIM map(2,2),co(2,2),iex(2),tex(2)
INPUT PROMPT"∃nをinputせよ n=":n
INPUT PROMPT"カオスの基準chao=":chao
FOR i=1 TO 2
   LET iex(i)=0
   LET tex(i)=0
NEXT i
READ co(1,1),co(1,2),co(2,1),co(2,2)
DATA 1,4,3,6
READ pai,ep1
DATA 3.141592653589793,.00001
REM ------
LINE INPUT PROMPT"テスト用(t) or 本番(h)?   t/h=":tq$
IF tq$="H" OR tq$="h" THEN GOTO 10
READ txl,txu,tyl,tyu,n1
DATA .3,.5,.6,.8,200
GOTO 20
10 LET txl=.1
   LET txu=.9
   LET tyl=.1
   LET tyu=.9
   LET n1=600
20 REM --------画像パラメタの設定--------------
   INPUT PROMPT"iteration の回数(80以上)=":tm
   PRINT "計算する四角形のinput"
   input PROMPT"(xl,xu)=":xl,xu
   INPUT PROMPT"(yl,yu)=":yl,yu
   PRINT"カラーの設定は以下の通り"
   PRINT"   O+(±√a)  共に収束； 黒 ｜O+(√a)収束,O+(-√a)non；赤"
   PRINT"O+(√a)non,O+(-√a)収束；緑 ｜O+(√a)non, O+(-√a)non；黄"
   INPUT PROMPT"カラーかえるか？(y/n)=":q$
   IF q$="N" OR q$="n" THEN GOTO 30
   print"黒はかえられません"
   INPUT PROMPT "co(1,2),co(2,1),co(2,2)をinputせよ":co(1,2),co(2,1),co(2,2)
   REM -----
30 LET ts=TIME
   LET wx=xu-xl
   LET wy=yu-yl
   LET twx=txu-txl
   LET twy=tyu-tyl
   LET dx=wx/n1
   LET dy=wy/n1
   LET a1=twx/wx
   LET a2=twy/wy
   LET b1=(xu*txl-xl*txu)/wx
   LET b2=(yu*tyl-yl*tyu)/wy
   PRINT "カラー   | O+(-√a);収束 | O+(-√a);non"
   PRINT "---------|---------------|---------------"
   IF q$="N" OR q$="n" THEN GOTO 40
   PRINT "√a;収束 |      黒       |    ";co(1,2)
   PRINT "√a;non  |      ";co(2,1);"      |    ";co(2,2)
   GOTO 45
40 PRINT "√a;収束 |      黒       |     赤"
   PRINT "√a;non  |      緑       |     黄"
45 print"-------------------------------------------"
   SET AREA COLOR 1
   PLOT AREA:txl,tyl;txu,tyl;txu,tyu;txl,tyu;txl,tyl
   PLOT TEXT, AT .02,.94:"Q2Mandel, aの範囲;"
   PLOT TEXT, AT .25,.94,USING"##.####### ##.#######":xl,xu
   PLOT TEXT, AT .55,.94,USING"##.####### ##.#######":yl,yu
   PLOT TEXT, AT .84,.94:"tm="   
   PLOT TEXT, AT .88,.94,USING"#####":tm
   REM
   REM -------------計算開始-----------------
   REM
   SET POINT STYLE 1 
   FOR i=0 TO n1
      LET ar=xl+i*dx
      FOR j=0 TO n1
         LET ai=yl+j*dy
         LET t=0
         LET r1=SQR(ar*ar+ai*ai)
         REM ------ 計算削除  --------------
         IF ar <0 THEN GOTO 70
         LET ia=INT(ar+1)
         IF ia=1 THEN GOTO 100
         IF ia=2 THEN GOTO 120
         IF ia=3 THEN GOTO 140
         IF ia=4 THEN GOTO 160
         GOTO 90
         REM
70       IF ABS(ai) >1.1 THEN GOTO 180
         LET zr=ar-.5
         LET zi=ai
         GOSUB 1000
         IF r > 1.54 THEN GOTO 180
         LET zr=ar+.5
         LET zi=ai
         GOSUB 1000
         IF r <0.45 THEN GOTO 180
         GOTO 90 
         REM
100       LET zr=ar-.5
          LET zi=ai
          GOSUB 1000
          IF r>1.4 THEN GOTO 180
          LET zr=ar*ar-ai*ai-3*ar+1
          LET zi=2*ar*ai-3*ai
          GOSUB 1000
          IF r<1 THEN GOTO 180
          GOTO 90
          REM
120       LET zr=ar-1
          LET zi=ai
          GOSUB 1000
          IF r<1 THEN GOTO 180
          LET zr=ar-.5
          LET zi=ai
          GOSUB 1000
          IF r>1.54 THEN GOTO 180
          GOTO 90
          REM
140       LET zr=ar-2.5
          LET zi=ai
          GOSUB 1000
          IF r>0.7 THEN GOTO 180
          LET zr=ar*ar-ai*ai-3*ar+1
          LET zi=2*ar*ai-3*ai
          GOSUB 1000
          IF r<1 THEN GOTO 180
          GOTO 90
          REM 
160       LET r5=ABS(ai)
          LET r6=-ar+3.5
          IF r5>r6 THEN GOTO 180
          LET zr=ar^3-2*ar*ar-3*ar-1+ai*ai*(2-3*ar)
          LET zi=ai*(3*ar*ar-ai*ai-4*ar-3)
          GOSUB 1000
          IF r<1 THEN GOTO 180
          GOTO 90
          REM
180       LET map(1,1)=map(1,1)+1
          GOTO 400
          REM ------ critical points の計算 ------
90       LET zr=ar
         LET zi=ai
         GOSUB 1000
         LET r=SQR(r)
         LET cx=r*COS(s/2)
         LET cy=r*SIN(s/2)
         REM------------------
         REM     loop1
         REM
         LET x=cx
         LET y=cy
200       LET d1=x*x-y*y-ar*x+ai*y+ar
          LET d2=2*x*y-ai*x-ar*y+ai
          LET d=d1*d1+d2*d2
          IF d<ep1 THEN LET iex(1)=iex(1)+1
          LET x1=(x*d1+y*d2)/d
          LET y1=(y*d1-x*d2)/d
          LET d3=SQR(x1*x1+y1*y1)
          LET t=t+1
          IF d3>chao THEN GOTO 205 ELSE GOTO 210
205       IF t>n THEN GOTO 250
          LET tex(1)=tex(1)+1
          LET x=x1
          LET y=y1
          GOTO 200
210       IF t<tm THEN GOTO 212 ELSE GOTO 215
212       LET x=x1
          LET y=y1
          GOTO 200
215       LET ii=1
          GOTO 290
250       LET ii=2
          REM
          REM        loop2
          REM
290       LET x=-cx
          LET y=-cy
          LET t=0
300       LET d1=x*x-y*y-ar*x+ai*y+ar
          LET d2=2*x*y-ai*x-ar*y+ai
          LET d=d1*d1+d2*d2
          IF d<ep1 THEN LET iex(2)=iex(2)+1
          LET x1=(x*d1+y*d2)/d
          LET y1=(y*d1-x*d2)/d
          LET d3=SQR(x1*x1+y1*y1)
          LET t=t+1
          IF d3>chao THEN GOTO 305 ELSE GOTO 310
305       IF t>n THEN GOTO 350
          LET tex(2)=tex(2)+1
          LET x=x1
          LET y=y1
          GOTO 300
310       IF t<tm THEN GOTO 312 ELSE GOTO 315
312       LET x=x1
          LET y=y1
          GOTO 300
315       LET jj=1
          GOTO 360
350       LET jj=2
          REM ----------  塗り絵 -------------
360       LET map(ii,jj)=map(ii,jj)+1
          IF ii=1 AND jj=1 THEN GOTO 400
          LET tx=a1*ar+b1
          LET ty=a2*ai+b2
          SET POINT COLOR co(ii,jj)
          PLOT POINTS :tx,ty
          REM --------------------------------
400    NEXT j
500 NEXT i
    REM
    REM
    REM
    LET tt=TIME-ts
    print"計算時間(sec=)";tt
    PRINT "----- カラーmap ------------"
    PRINT"(1,1), (1,2)=";map(1,1),map(1,2)
    PRINT"(2,1), (2,2)=";map(2,1),map(2,2)
    print"------------------------------"
    PRINT"反復∃n以下,絶対値chao以上の回数=";tex(1),tex(2)
    PRINT"分母ほとんど0の回数=";iex(1),iex(2)
    REM
    PLOT TEXT, AT .02,.05:"カラーmap="
    PLOT TEXT, AT .16,.05,USING"###### ######  ###### ######":map(1,1),map(1,2),map(2,1),map(2,2)
    PLOT TEXT , AT .6,.05:"∃n, chao="
    PLOT TEXT, AT .73,.05,USING"#### ###.##":n,chao
    LINE INPUT PROMPT"ファイル名など書くか？ y/n ":fa$
    IF fa$="N" OR fa$="n" THEN GOTO 1100
    PRINT"ファイル名などinputせよ，(例) Q2m-01,Q2m-02..."
    LINE INPUT PROMPT"ファイル名など（全角20文字位）：":na$
    PLOT TEXT, AT .05,.01:"file名など:"
    PLOT TEXT, AT .2,.01:na$
    print"***最後に絵のセーブ忘れずに***"
    GOTO 1100
    REM ==============================================================
1000 REM--- サブルーチン polor; zr,zi をinputして r,θ(s) を返す
     REM   円周率πの値 pai=3.141592... はメインプリグラムに必要
     REM ==============================================================
     LET zero=1.e-7
     LET r=SQR(zr*zr+zi*zi)
     IF ABS(zr)>= zero THEN GOTO 1010
     IF zi>0 THEN 1020 ELSE 1030
1020 LET s=pai/2
     RETURN
1030 IF zi<0 THEN 1040 ELSE 1050
1040 LET s=3*pai/2
     RETURN
1050 LET s=0
     RETURN
     REM ----------- 
1010 LET AT=ATN(zi/zr)
     LET s=AT+pai
     IF zr>0 THEN LET s=AT
     RETURN
     REM =====================================
1100 END